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ABSTRACT 

Fingering convection (otherwise known as thermohaline convection) is an 
instability that occnrs in stellar radiative interiors in the presence of nnstable 
compositional gradients. Nnmerical simnlations have been nsed in order to esti¬ 
mate the efficiency of mixing indnced by this instability. However, fnlly three- 
dimensional (3D) compntations in the parameter regime appropriate for stellar 
astrophysics (i.e. low Prandtl nnmber) are prohibitively expensive. This raises 
the qnestion of whether two-dimensional (2D) simnlations conld be nsed instead 
to achieve the same goals. In this work, we address this issne by comparing 
the ontcome of 2D and 3D simnlations of hngering convection at low Prandtl 
nnmber. We hnd that 2D simnlations are never appropriate. However, we also 
hnd that the reqnired 3D compntational domain does not have to be very wide: 
the third dimension need only contain a minimnm of two wavelengths of the 
fastest-growing linearly nnstable mode to captnre the essentially 3D dynamics of 
small-scale hngering. Narrow domains, however, shonld still be nsed with cantion 
since they conld limit the snbseqnent development of any large-scale dynamics 
typically associated with hngering convection. 

Subject headings: hydrodynamics - instabilities - stars : interiors - stars : evo- 
Intion 


Introduction 


Fingering convection (also called thermohaline convection) was hrst discnssed in the 
astrophysical context by Ulrich] (1972) (see also Ulrich||1971 ). The basic hngering instability, 
which occnrs in systems that are compositionally nnstably stratihed and thermally stably 
stratihed, is common in stars whose snrface is pollnted by high mean molecnlar weight 
material. This can happen throngh the accretion of material from infalling planets or debris 
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disks ( 

Vauclair 

2004 Garaud 20111 

Deal et al. 

2013 

) or from more evolved companion stars 

(Stanclihe et al 

2007 Angelou et al. 

2012 

). Fingering convection can also result from the 


radiative levitation of heavy elements such as iron from deeper regions upwards in relatively 
high-mass stars (|Theado et ah 2009; Zemskova et ah 2014), or from off-center nuclear burning 


in post main sequence stars (Charbonnel & Zahn 2007 Denissenkov 2010 Denissenkov & 
Merryheld 2011[ ). 


In the typical conditions encountered in stellar interiors, instability to hngering convec¬ 
tion only depends on the value of the density ratio Rq, dehned as 


Rr, — 


/ dTp _ \ 

v dr _ dr ) 

^*0 


( 1 ) 


where To(r) is the local temperature prohle in the stellar region considered, T^^{r) is the 
temperature prohle that region would have were it to be adiabatically stratihed, /io(?") is the 
local mean molecular weight prohle, and a and /3 are derivatives of the equation of state 
dehned as 

The density ratio is therefore the ratio of the density gradient due to the stabilizing thermal 
stratihcation to the density gradient due to the destabilizing compositional stratihcation. 
Note that Rq is also commonly written as 
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where V, Vad and have their usual astrophysical dehnitions, and 
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when 


As shown by Baines & Gill (1969), linear instability to hngering convection only occurs 

( 5 ) 


l<Ro< i?crit = r 1 


Kt 


K,. 


where kt and are the thermal and compositional dihusivities respectively. The lower limit 
{Rq = 1) corresponds to the onset of standard overturning convection, while the upper limit 
{Rq = Rcrit) corresponds to the critical point of marginal stability to hngering convection. 
As shown by Garaud et ah (2015), r varies between 10“® and 10“® in the interiors of main 
sequence stars, and increases only up to about 10“^ in degenerate regions of more evolved 
stars (e.g. red giants or white dwarfs). This shows that Rcrit is always much greater than one. 
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so that stellar interiors can very easily be destabilized: even a tiny adverse compositional 
gradient can trigger fingering instabilities. As a result, fingering convection is very common, 
and should be taken into account when considering mixing in stellar evolution models. 

An important difference in the behavior of fingering instabilities in astrophysical and 
geophysical systems comes from the value of the Prandtl number (Pr = zz/k-t, where v is the 
kinematic viscosity): this parameter is usually larger than one in most geophysical fluids, 
but is asymptotically small (typically, at most one order of magnitude larger than r; see 


Garaud et ah 2015) in stellar interiors. As a result, fingering structures which are typically 


fairly laminar in the geophysical context are instead very turbulent in the astrophysical one. 
This complicates the study of the saturation of the instability, which is necessary to estimate 
the rate of heat and compositional transport induced by fingering convection. 

In the past few years, thanks to the development of fast parallel algorithms and advances 
in supercomputing, it has become possible to follow the nonlinear development of fingering 
instabilities in three dimensional (3D) simulations at relatively low diffusivity ratio and low 
Prandtl numbers, albeit not yet at actual stellar values of these parameters. The first of 


such studies were by Denissenkov (2010) using 2D simulations and by Traxler et ah (2011) 


and Brown et ah (2013) in 3D. The 3D simulations covered a wide range of parameter space 


with Pr and r varying from 0.3 down to 0.01, and the density ratio was varied across the 
instability range (from 1 to r“^) 


Using them. Brown et al. (2013) were able to propose 
an analytical model for turbulent transport by fingering convection at low Pr and r. This 
model has only one free parameter, which can be calibrated using the 3D simulations. Once 


calibrated, the Brown et al. (2013) model correctly predicts the transport rate for both heat 


and composition for all available simulations in which r < Pr within a factor of 2 at worst, 
often much better. 

Whether this model remains as accurate for Pr and r lower than 10“^ still remains to 
be determined. However, cubic-domain 3D simulations at very low Pr and r rapidly become 
computationally prohibitive since the resolution must be increased as Pr and r decrease, in 
order to fully resolve the various boundary layers that form at the edge of the fingers. In 


addition, and as discussed by Traxler et al. (2011b) and Garaud et al. (2015), secondary 


large-scale instabilities can develop spontaneously from fingering convection, and the latter 
substantially affect the turbulent transport properties of the system. To study them requires 
very large computational domains, containing at the very least 20 to 30 wavelengths of the 
fastest-growing fingering mode in both horizontal and vertical directions, and often much 
more than that. Again, using such large domains in 3D, together with low Pr and r, is 
currently computationally prohibitive. 


For both reasons, it is very tempting to go back to 2D simulations either for very 


























4 


low diffusivity studies, or for very large domain studies. In fact, 2D fingering convection 
simulations used to be common in the oceanographic literature only 20 years ago (e.g. Shen 


1995 

Stern et ah 

2001 

), and are still often in use today ( 

Simeonov & Stern 

2007; Radko 

2008 

Sreenivas et al. 

2009 

Singh & Srinivasan 

2014 

) in that context. 

Stern et al. 

(2001) 


showed by comparing 2D and 3D simulations at Pr = 7 that the typical turbulent fluxes 
estimated from 2D simulations do not vary from those obtained in 3D simulations by more 


than a factor of a few. Denissenkov & Merryheld (2011) claim to have run both 2D and 3D 
simulations of hngering convection for Pr and r both 0(10“®), and to have found that the 
effective turbulent diffusivities of composition are very similar in the two cases. 


However, the work of Radko (2010) casts strong doubts on the claims of Denissenkov & 


Merryheld (2011). Indeed, he showed both analytically and numerically that the saturation 


of the hngering instability at low Pr in 2D can be attributed to the development of strong 
shear layers that destroy the hngers and completely change the overall behavior of the induced 


turbulence. A similar process was not seen in the 3D simulations of Radko & Smith (2012) 


and Brown et ah (2013). This raises a number of practical questions: (1) Is shear indeed 


present in low Prandtl number hngering (as claimed by Radko 2010), or not (as implied 


by Denissenkov & Merryheld 2011)? (2) More generally, under which circumstances, if any, 
can 2D simulations be used to study hngering convection at low Prandtl number? (3) If 
2D simulations are not appropriate, then how thick does a domain have to be in order to 
approximate 3D results appropriately? In this work, we answer all three questions. 


Section outlines our model equations, boundary conditions, and domain geometry. 
In Section we hrst compare 2D and 3D results in high density ratio simulations. In 
Section we present a similar comparison but in low density ratio simulations. As we shall 
demonstrate, in both cases 2D and 3D results are fundamentally diherent, but for diherent 
reasons. Section discusses our hndings and proposes acceptable compromises in terms of 
running 3D hngering simulations as cheaply as possible. 


2. Model 


In what follows, we study the dynamics of homogeneous hngering convection in the 


Boussinesq approximation, as introduced for the purpose of numerical simulations by |Shen 
(1995) (see also Radko 2003). In this formalism, the background temperature and compo¬ 


sition prohles are assumed to be linear, with locally constant gradients dTo/dr and dfio/dr 
(see Section]^, while all perturbations to that background (including the temperature T, 
the pressure p, the velocity held u = (u, u, w) and the mean molecular weight /i) are assumed 
to be triply periodic in the domain considered. 
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The non-dimensional governing equations are given by: 


1 f \ 

+ u ■ VuJ = -Vp + (T - /i)e, + V'u, 

dT 

— + u ■ VT + u; = V^T, 

ot 

d^i 1 2 

^ + u ■ V/i + —w = 
ut Rq 

V-u = 0 


( 6 ) 


where Pr, r, and Rq were introduced in Section To arrive at these equations, we have 
used the following units: 


1/4 


[i] = d = 




Q/Q I *^^0 dTa.d I 

^ \ dr dr 


as the unit distance, 


d? kt 

[t] = —, [m] = — as the unit time and velocity respectively. 


[r] = 


dTn dTf 


dr 


ad 


dr 
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d as the unit temperature, [/i] = — [T] as the unit composition) 


where g is the local gravity, and all other quantities were dehned in Section 

In what follows we consider two possible types of domains: 2D domains in the (x, z) 
plane of size (250d x 250d), and 3D domains of size (250d x Ly x 250d), where Ly shall be 
varied. We use a fairly large size in the {x, z) plane to allow for the development of large-scale 
structures in the simulations, should they naturally occur. We set the parameters Pr and r 
to be both equal to 0.03, well below unity. The density ratio, on the other hand, is varied to 
detect and study the existence of different dynamical regimes across the hngering instability 
range as appropriate. The code used to solve the set of equations ([^ under triply-periodic 
boundary condition is the pseudo-spectral code described in [Traxler et ah (2011b). In all 
cases, we use a resolution of 256 spectral modes per 250d (corresponding to an effective 
resolution of about 3 grid points per d). The same resolution is used in all directions. 


3. Comparison of 2D and 3D runs at high density ratio 


We begin our investigation by considering a “high” density ratio case, with Rq = 5. 
While still signihcantly below i?crit = 33.3 (at these parameters), this value is already large 
enough to be in the sheared regime described by Radkoj (2010) (see Section [^, but small 
enough to ensure that resolving the buoyancy frequency (which is equal to iyPr(i?o ~ 1) iii 
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these units) remains computationally manageable. Indeed, gravity waves are present in this 
system, and the time step must be small enough to fully resolve their oscillation period. 


3.1. Results of 2D simulations 


Figure [T] shows the evolution of the mean molecular weight perturbation fi and of the 
horizontal component of the velocity held u in the 2D run, from very early times to much 
later times. We see, successively, the development of the basic hngering instability (top), its 
saturation (middle), and hnally, the emergence of strong alternating shear layers and their 
effect on the hngers (bottom). As seen in Figure]^ the shear layers evolve slowly with time: 
weak layers appear to be pushed closer to one another, and eventually merge into stronger 
ones. As discussed in Section 1, the spontaneous emergence of large-scale shear layers in 
2D hngering convection at low Prandtl number had already been discovered and studied by 


Radko (2010). We therefore conhrm his hndings here. 


Interestingly, we further hnd that the nonlinear interaction of the hngers and the shear 
drive what appears to be relaxation-oscillations. The latter are most noticeable in Figure 
when viewed in terms of the total compositional hux {wfi) and of the r.m.s. horizontal 
velocity Urms = (where the brackets (■) denotes a volume average over the entire do¬ 

main). The oscillations appear around t = 5000, then grow in amplitude as the various shear 
layers merge. By t = 15000, the oscillation pattern becomes very clear: efficient hngering 
(characterized by a large vertical hux) drives horizontal shear (noticeable in the increase of 
Mrms), which hrst distorts the hngers then eventually completely suppresses hngering con¬ 
vection (hence the drop in the turbulent hux). When this happens, the mechanism that 
drives the shear dies out, and the shear gradually disappears. The hngers are eventually 
allowed to grow again, and the cycle repeats. Note that the amplitude and period of the 
cycle appear to depend on the number of shear layers present. As the simulation progresses 
and the number of layers decreases through mergers, the cycle lengthens and its amplitude 
grows, until a quasi-periodic state is reached. 


3.2. Results of 3D simulations of various Ly 

Even though these fascinating dynamics appear in 2D simulations of hngering convec¬ 
tion, it is easy to show that the latter are spurious and do not exist in computational domains 
that are “sufficiently 3D”. Figure shows the turbulent compositional hux (tn/i) and the 
r.m.s. horizontal velocity Urms as a function of time for simulations that are perfomed at 
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Fig. 1.— Snapshots of the compositional pertnrbations (left) and of the horizontal velocity 
in the x direction (right) at three times (see the evolntion also shown in Figure]^: near the 
onset of the hngering instability {t = 280, top), early after its satnration {t = 530, middle) 
and at later time once the shear layers have begun to develop {t = 2500 bottom). Note how 
the horizontal shear distorts the hngers. 























- 8 - 


1.5 


1 

0.5 

0 


-0.5 


-1 

-1.5 

0 10000 20000 30000 40000 50000 

time 

Fig. 2.— Horizontally-averaged horizontal velocity profile in the 2D simulation of hngering 
at i?o = 5 as a function of time. The layers rapidly merge down to four, after which a 
quasi-steady periodic state appears to be reached. The oscillations in the intensity of the 
shear are clearly visible. 




Fig. 3.— Compositional flux and r.m.s. horizontal velocity in the 2D simulation of hngering 
at Ro = 5, illustrating the mechanism responsible for the relaxation oscillations. 
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exactly the same parameters and have Rq = 5, as in Section 3.1, but now in 3D with differ¬ 
ent domain thicknesses Ly. The early stages of the 2D simulation discussed in the previous 
section are included for comparison. 


It is quite clear that the simulation for which Ly = Ad still behaves as if it were 2D, that 
is, with a turbulent flux that rapidly decreases with time after the initial saturation, and a 
corresponding r.m.s. horizontal velocity that increases with the generation of shear layers. 
By contrast, simulations with Ly = 15d and larger all behave in quantitatively similar ways, 
achieving the same well-defined quasi-steady state post-saturation, with substantially higher 
flux and lower horizontal velocities than in 2D and no obvious quasi-periodic oscillations. 
The run with Ly = 8d appears to have the same character as the wider 3D domains but 
somewhat underestimates the flux. 


The significance of the cutoff-sizes Ly = 8d and Ly = 15d becomes clearer if one notes 
that, at the parameter values selected, the typical wavelength of the fastest-growing mode of 
fingering convection is about 8.5d (one wavelength containing one up-going and one down¬ 
going finger). We conclude from this experiment that at low Prandtl number and high 
density ratio: (1) a domain has to contain at least 1 wavelength of the fastest-growing 
hngering mode in the third dimension to exhibit dynamics that are fully 3D instead of the 
spurious shear layers and associated relaxation oscillations that are observed to exist in 2D; 
(2) a domain has to contain at least 2 wavelengths of the fastest-growing mode to yield 
quantitatively accurate results for the turbulent fluxes and turbulent intensity of fingering 
convection. These results are generic: similar quantitative findings have been obtained for 
other values of Pr, r <C 1 and for sufficiently large Rq. 


3.3. Can the artificial shear be prevented? 


The results of Figure strongly suggest that the substantial horizontal shear and the 
associated quenching of the vertical fluxes observed in the 2D case are artificial. However, 
given the vast savings in computational time between 2D and 3D simulations, it is tempting 
to wonder if one could still recover acceptable 2D solutions (i.e. solutions whose dynamical 
behavior and turbulent transport rate are comparable to the 3D solutions) by suppressing 
the shear. Figures]^ and [^demonstrate that this does not work. Since the code we are using 
is spectral, it is easy to zero out the mean component of the horizontal flow at every time 
step, effectively suppressing the shear whilst allowing all other modes to evolve normally. We 
performed this experiment on the 2D case at i?o = 5 discussed in Section 3.1[ We compare 
below these new results to those originally obtained without shear-suppression, and to the 


narrow 3D case with Ly = 15d discussed in Section 3.2 that appeared representative of the 
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true solution. Finally, in order to ensure that shear-suppression does not artificially affect 
the dynamics of fingering convection in 3D, we ran a hnal experiment suppressing the shear 
in the same manner on another wider 3D case (with domain size 250d x 50d x 250d) at the 
same parameters. 

Figure [^shows the thermal flux (wT) for each of these cases. As seen before, the narrow 
3D case settles down to a much higher flux than the 2D case, which exhibits relaxation 
oscillations. The wider 3D case with the shear suppressed has the same flux as the narrower 
3D one, showing that shear does not play a signihcant role in the 3D simulations. The 2D 
case with the shear suppressed, however, does something altogether different from all the 
other cases. Although it initially settles into a mode with roughly the right mean flux, the 
fluctuations about the mean are much larger than in the 3D cases. Furthermore, at around 
t = 40000, the 2D shear-suppressed case jumps to a solution with a much higher flux. 


This rather peculiar behavior is best understood by looking at snapshots of the simula¬ 
tions. Figure]^ which shows slices of the horizontal velocity field, demonstrates the physical 
differences between these various cases. In 3D, the solution consists of a homogeneous field 
of small-scale hngers, whether the shear is removed or not. In the original 2D case, the pres¬ 
ence of strong shear discussed in Section 3T is obvious. In 2D with the shear suppressed, it 
appears that the solution attempts to create a velocity pattern that is as close to the regular 
2D case (with strong mean shear) as possible, given the model constraints. After a short 
transient phase, strong bands of flow appear. They are not exactly horizontal, as that is 
disallowed, but are slightly angled upwards instead. To begin with, the specific angles are 
somewhat random over space. The corresponding vertical flux is significantly larger than 
in the regular 2D case, and somewhat comparable to the one observed in 3D. Ultimately 
though, the flow self-organizes at a single well-defined angle, and acts as an extremely effi¬ 
cient periodic “escalator” that significantly enhances the vertical transport. Despite being 
quite interesting, however, these flows remain clearly artihcial compared to the 3D cases. 
We are therefore forced to conclude that there is no easy £x to the problem, and that 3D 
simulations are the only way to get reliable results. 


4. Comparison of 2D and 3D runs at low density ratio 

We now turn to low density ratio systems, and illustrate our findings with a case which 
has Rq = 1/0.9 ~ 1.11, and Pr = r = 0.03 as above. Figure shows snapshots of a 2D 
simulation, and of a 3D simulation done in the domain of thickness Ly = 125d, both taken 
around t = 500. Again, we see a notable difference between the 2D and the 3D case, but 
the nature of this difference is clearly not the same as the one discussed in Section In the 
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2D case, we see large plumes that are reminiscent of oversized turbulent fingers. It is worth 
remembering that the typical basic fingering mode width is about 1/30 of the domain size, 
so the structures seen here (which have a size of about 1/3 of the domain) are much larger 
than basic fingers. The 3D case, by contrast, exhibits very different dynamics. We note a 
clear separation of scales between small-scale fingers (which have more or less the same size 
as the basic instability) and domain-scale gravity waves that modulate the fingering field. 


The emergence of these large-scale gravity waves is well-understood, and can be at¬ 
tributed to the collective instability hrst discussed by Stern (1969) in the oceanographic 


context, later found in simulations of low-Prandtl number hngering convection by Brown 


et al. (2013), and systematically studied by Garaud et ah (2015). The collective instability 


is a mean-held instability, excited through a positive feedback loop between the large-scale 
gradients of temperature and composition, and the respective turbulent huxes induced by 
these large-scale gradient^ These large-scale gravity waves do not appear to be excited 
in 2D, which is somewhat surprising since the collective instability theory does not a priori 
need to be 3D to operate. 


A more quantitative way of seeing the scale separation inherent to the collective insta¬ 
bility in 3D, and its absence in 2D, is by inspection of the kinetic energy spectrum. This 
is shown in Figure as a function of the vertical wavenumber, for the same simulations 
and at the same times as the ones shown in Figure The spectrum of the 2D run clearly 
peaks around k = O.OSd”^, which corresponds to a typical wavelength of about 80d, or in 
other words, about 10 times the wavelength of the fastest-growing hngering mode. This 
is more-or-less the scale of the plumes observed in the snapshot. Between that value of 
k and the energy injection scale {k ~ 1) we observe that the kinetic energy spectrum is 
close to a power law with index —2, which may indicate the presence of an inverse energy 
cascade. The 3D run, by contrast, has a kinetic energy spectrum that peaks at the lowest 
possible wavenumber k = 0.025d“^, which corresponds to a wavelength commensurate with 
the domain size. This is indeed the scale of the collective instability mode observed in the 
snapshot. The energy spectrum drops sharply between k = 0.025d“^ and k = 0.1d“^, then 
more gently between k ~ O.ld”^ and k ^ d (the energy injection scale). Beyond k = d, both 
2D and 3D simulations have the same energy spectrum which is dominated by small-scale 
diffusive processes, showing that they are basically the same in 2D and 3D. 

Finally, Figure shows the turbulent compositional flux for simulations where the thick- 


^As discussed by |Radko| ( |2013[ ), the collective instability can be interpreted as the turbulent analog of 
Oscillatory Double-Diffusive Convection, a linear instability that normally takes place in semiconvective 
regions ( Kato|[l9^ ). 
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ness of the domain Ly is varied from 0 to 125(i as we did in Figure for the high density 
ratio case. This time, we see a clear difference in behavior between all the 3D simulations 
and the 2D simulation. The 2D case shows a continued growth of the efficiency of turbulent 
transport that only saturates quite late (after t = 500) at a fairly high mean value. All the 
3D simulations, by contrast, show an early saturation of the primary fingering instability 
around t = 150, then a secondary growth phase associated with the growth of the collective 
instability, which later also saturates but at a much lower mean value than the 2D case. 
By contrast with the high density ratio case described in the previous Section, even the 
Ly = Ad simulation seems to be qualitatively more compatible with the full 3D case than the 
2D case. However, in order to get quantitatively accurate flux measurements just after the 
saturation of the hngering instability but prior to the onset of the collective instability (i.e. 
roughly between t = 150 and t = 250), it is necessary to choose Ly > 15d, just as in the high 
density ratio case. In other words, simulations in domains that are at least 2 wavelengths of 
the fastest-growing mode (which is also about 8.5d at this value of Rq) are necessary to get 
adequate estimates of the fingering fluxes at low Prandtl number. 


Of course, even if the fingering dynamics are correctly accounted for, a very narrow 
domain may no longer be adequate once larger-scale structures appear. In the Ly = 15d 
domain, the large-scale gravity waves excited by the collective instability are artihcially 
forced to be 2D, for instance. How this affects their saturation amplitude, and their overall 
behavior, remains to be determined. Similarly, in the work of Zemskova et ah (2014), the thin 
domain used severely restricts the dynamics of the convective layer that eventually forms, 
to the extent that the results in the convective phase are not necessarily reliabl^ In other 
words, while a thin domain is a good compromise to study homogeneous fingering convection, 
a full 3D domain with aspect ratio of order unity may remain necessary to properly model 
any large-scale dynamics that naturally arise. 


5. Discussion and prospect 

In this work we have compared 2D and 3D simulations of fingering convection at low 
Prandtl number to determine under which circumstance 2D models can be used to approx¬ 
imate 3D systems. We have found that 2D simulations always suffer from some kind of 
pathology, but that pathology is different for different density ratios. 

For more strongly stratified systems (higher density ratio), 2D simulations produce arti- 


^see, for instance, the strong vortices that are clearly visible in their Figure 10, a classic sign of quasi-2D 
dynamics that doesn’t persist in 3D 
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ficial shear layers that strongly suppress the efficiency of hngering convection, and nonlinearly 
interact with the latter to produce relaxation oscillations. The emergence of shear at high 
density ratios in 2D does not come as a surprise, as it was predicted and demonstrated to 
occur by Radko (2010). Moreover, a similar phenomenon occurs in 2D Rayleigh-Benard con¬ 
vection in vertically-bounded but horizontally-periodic models (Goluskin et ah 2014). The 
3D simulations, by contrast, do not show any signihcant shear. The turbulent fluxes in 3D 
are signihcantly higher than in 2D, and are more-or-less steady once the system has reached 
saturation. 


For more weakly stratihed systems (lower density ratio), 2D simulations appear to ex¬ 
hibit an inverse energy cascade which results in the progressive coalescence of hngers into 
larger and larger ones, a phenomenon that is not seen in 3D. As a result, the fluxes at satu¬ 
ration are signihcantly higher in 2D than in 3D. Furthermore, this inverse cascade seems to 
prevent the development of the collective instability, which in 3D drives large-scale gravity 


waves (Garaud et al. 2015). 


In short, 2D simulations always fail to reproduce the basic dynamics of hngering con¬ 
vection at low Prandtl number and should therefore never be used in this context. The fact 


that they appear to be adequate for high-Prandtl number studies (Stern et al. 2001) is, in 


the light of our work, somewhat surprising. However, this could be due to the fact that the 
relative importance of the nonlinear terms in the momentum equation is inversely propor¬ 
tional to the Prandtl number (see equation]^, and that in both cases, the ohending artihcial 
behavior discovered in 2D is an inherently nonlinear one. A preliminary study reveals that 
the transition from dynamically-correct to pathologically-sheared 2D hngering convection 
occurs for Prandtl numbers around 0.5; in other words, for Pr > 0.5, 2D simulations provide 
qualitatively reasonable results, but for Pr < 0.5, 3D simulations are necessary. 


Our hndings cast doubt on the 2D hngering hux estimates reported by |Denissenkov 
(2010) and Denissenkov & Merryheld ( 2011[ ). In fact, the simulation snapshots in Figure 
4 of Denissenkov (2010) are quite similar to our 2D snapshots at t = 280 in Figure 
and clearly show the early stages of the shear development with strongly distorted hngers. 


Meanwhile, the claim made by Denissenkov & Merryheld (2011) that 2D and 3D huxes at 


low Prandtl number are very similar could be attributed to the fact that they only integrated 
their simulations for very short times: we hnd here that 2D and 3D huxes are closer to one 
another before the shear has gained signihcant amplitude. 


We have also found, however, that a fully 3D domain (with equal size in all dimen¬ 
sions) is not necessary to model the correct behavior: a minimum domain width of about 
2 wavelengths of the fastest-growing hngering mode is sufficient to eliminate the unwanted 
behavior, and to yield estimates of the turbulent huxes that are quantitatively consistent 
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with those obtained in nearly cubic domains. Narrow domains, however, should still be used 
with caution since they could limit the subsequent development of any large-scale dynamics 
typically associated with fingering convection (e.g. layer formation and the excitation of 
large-scale gravity waves). 


Interestingly, while there is a fundamental difference between 2D and 3D fingering sim¬ 
ulations, Moll et ah (2015) found that this is not the case for simulations in the oscillatory 
double-diffusive regime (i.e. in the semiconvective regime) at similarly low Prandtl num¬ 
bers. There, 2D and 3D simulations behave in qualitatively similar ways, and the fluxes 
extracted in 2D and 3D are within an order of magnitude of each other for all parameters 
surveyed. This striking difference between the hngering regime and the oscillatory regime 
raises a number of questions. Why are they so different, given that the only difference in the 
governing equations is the sign of the background temperature and compositional gradients? 
Is the “two-wavelength” guideline we propose for the selection of a minimal 3D domain size 
in hngering convection exportable to other instabilities? Generally speaking, can one pre¬ 
dict ahead of time and without running simulations whether a given set of equations and 
boundary conditions, solved for a given set of parameters in 2D, will be a good approxima¬ 
tion to the full 3D dynamics? These are clearly difficult mathematical questions, but their 
answers, if they exist, could provide formal guidelines for the general reliability of 2D and 
narrow-domain 3D simulations. 
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Fig. 4.— Compositional flux (left) and r.m.s. horizontal velocity (right) for hngering simu¬ 
lations at Ro = 5, of varying domain thicknesses (see legend for detail). Runs with very thin 
domain behave as if they were 2D, while thick enough domains behave as in a 3D manner 
that rapidly becomes independent of the domain size. 
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Fig. 5.— Comparison of time evolution of thermal fluxes for 2D and 3D simulations at 
Ro = 5 with and without mean shear flows suppressed. Fluxes are shown for a purely 
2D simulation, a purely 2D simulation with shear suppressed, a narrow 3D simulation that 
contains 2 hnger wavelengths in the narrow direction, and a wider 3D simulation with the 
shear suppressed. Clearly, whilst both 3D simulations agree well, both 2D simulations do 
not reflect the same behavior at all. 
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Fig. 6.— Slices in the (x, z) plane of the horizontal velocity for the simulations shown in 
Fig. 1^ at late times in each simulation. In 2D, strong shear flows are seen, and in 2D with 
mean shear suppressed, an “elevator” solution is observed that is very close to the sheared 
2D case but which is periodic. None of this behavior is seen in 3D, where the suppression of 
shear makes no difference since signihcant shear is not present. 



Fig. 7.— Snapshots of the compositional perturbation for a simulation with Rq = 1.11, in 
2D (left) and in a thick 3D domain (right). Note the presence of large-scale gravity waves 
in the 3D case, absent in the 2D case. Also note the signihcant difference in the maximum 
amplitude of the compositional perturbations in each case. 
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Fig. 8.— Kinetic energy spectrum as a function of the vertical wavenumber k, for the two 
simulations shown in Figure The linear instability injection scale is estimated from the 
horizontal wavenumber I of the fastest-growing mode, and assuming that k ^ 1. 
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Fig. 9.— Compositional flux for Rq = 1.11, for various domain thicknesses (see legend for 
detail). The 2D and 3D runs behave very differently. 














